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The first and second order hyperpolarizabilities have been extensively studied to identify universal 
properties near the fundamental limit. Here, we employ the Monte Carlo method to study the 
fundamental limit of the second hyperpolarizability. As it was found for the first hyperpolarizability, 
the largest values of the second hyperpolarizability approaches the calculated fundamental limit. 
The character of transition moments and energies of the energy eigenstates are investigated near 
the second hyperpolarizability's upper bounds using the missing state analysis, which assesses the 
role of each pair of states in their contribution. In agreement with the three-level ansatz, our 
results indicate that only three states (ground and two excited states) dominate when the second 
hyperpolarizability is near the limit. 
OCIS Codes: 190.0190, 020.0020, 020.4900 



I. INTRODUCTION 

In 2000, the fundamental limit of the off-resonant elec- 
tronic first hyperpolarizability and second hyperpolariz- 
ability were calculated in an attempt to understand con- 
straints imposed by quantum mechanics [H-Q ■ (Note that 
for the rest of this paper, we say hyperpolarizability when 
referring to the first hyperpolarizability.) Using Thomas- 
Kuhn sum rules and assuming that when the hyperpolar- 
izability ft is optimized only three states (including the 
ground state) contribute, an upper limit of ft was found 
to be given by, 

where — e is charge of the electron, Ti is Planck's con- 
stant, N denotes the number of electrons in the atom or 
molecule and £m is the energy difference between the 
first excited state and ground state. In conjugated sys- 
tems, it is known that only the 7r-electrons contribute 
to the nonlinear-optical response, so the effective num- 
ber of electrons, N = N e ff is used instead of the total 
number. 0, [E[ 

A comparison the largest experimentally measured sec- 
ond order molecular susceptibilities with the fundamental 
limit reveals a large gap between the two. Measurements 
have never crossed the limit and are typically well below 
it. Until about 2005, all molecules ever measured mys- 
teriously fell a factor of about 30 below the limit, @, 0] 
until a new class of molecules were reported by Kang and 
coworkers. [H Q 

In early studies, several experimental tools such as 
linear spectroscopy, Raman spectroscopy, hyper- Rayleigh 
scattering and Stark spectroscopy were used to inves- 
tigate the possible reasons for the gap between theory 
and experiment. There it was explained that the gap 
is mainly due to the unfavorable arrangement of excited 
state energies. 

In an effort to understand how one might make a real 
quantum system that has a hyperpolarizability near the 



limit, Zhou and coworkers applied a numerical approach 
to determine the shape of the potential energy function 
that yields the largest value of the hyperpolarizability [111 
[l2j . The potentials were chosen among a wide range of 
possibilities including polynomials, power laws and piece- 
wise continuous functions. The best hyperpolarizabili- 
ties found were given by ft sa 0-7ft max . Thus, a new but 
smaller gap appeared between numerical simulations and 
the fundamental limit. These studies suggested that the 
concept of conjugation of modulation could be used as 
a new design paradigm for making better molecules that 
was later experimentally verified with one-dimensional 
molecules that broke above the apparent limit. fl3l [l4| 

Those numerical optimization studies showed the first 
hint of the appearance of universal properties of quan- 
tum systems that approach the fundamental limit. In all 
cases, independent of the starting potential, a large range 
of optimized systems were found to have (1) ft sa 0.7ft max ; 

(2) they are all well approximated by a three-level model; 

(3) the energy ratio between the first two excited states is 
always near 1/2; and, (4) the the normalized transition 
moment to the first excited state is about 0.75. [HI [l2[ 
A broader set of studies found similar universal proper- 
ties, for example, when considering the effects of geom- 
etry (i.e. the optimal arrangement of point nuclei) fl5| 
and the effects of an externally applied electromagnetic 
field. 

Since the 30% gap may be an indication of a flaw in 
the calculation of the fundamental limits, Kuzyk and 
Kuzyk (KK) focused on investigating the validity of the 
three level ansatz, the only assumption underlying limit 
theory. [Tt} The three level ansatz states that when the 
hyperpolarizability of a quantum system is at the fun- 
damental limit, only three states contribute; Q that is, 
the contributions from all other states is negligible. Note 
that the ansatz does not claim the converse. When three 
states are involved, the hyperpolarizability is not neces- 
sarily at the limit. KK used Monte Carlo simulations to 
numerically generate a set of ft values to investigate if 
any of the data attains the predicted fundamental limit. 
They found a continuous distribution of values up to the 
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fundamental limit. These studies verified the validity of 
the theory of limits and therefore indirectly verified the 
three level ansatz. 

More interestingly, these results imply that a set of 
wavefunctions exist that are consistent with the sum 
rules, yet may not be derivable from typical Hamiltoni- 
ans; that is, the sum rules may permit state vectors that 
are not derivable from the Schrodinger equation. This 
statement may seem contradictory since sum rules are 
directly derivable from the Schrodinger equation; but, 
the sum rules apply under broader conditions, which we 
can understand as follows. 

The sum rules are derived using the fact that 
[x, [x, H]] = Ti 2 /2m for a typical Hamiltonian. How- 
ever, it may be possible that [x, [x, H + A]] — h 2 /2m 
for an operator A, where A can take a form that is 
never observed in molecular systems. As an example, 
it is straightforward to show that the unusual operator 
A = /(x, y, z)p x PyPz obeys the commutation relationship 
yet is of a form that is to our best knowledge never found 
in nature. Other such operators can be invented. Thus, 
while all standard Hamiltonians obey the sum rules, more 
exotic Hamiltonians can also obey the sum rules. 

Sum rules have been used to determine un- 
known transition moments when fitting the dispersion 
data of the second harmonic hvperpolarizabilitv. (l8l| 



and for understanding length-scaling of the second 
hyperpolarizability. [l9j Fundamental limits have also 
been used to investigate the essential ingredients re- 
quired of small molecules to approach the fundamental 
limit. (20l. l2~l| However, no numerical optimization stud- 
ies have targeted the second hyperpolarizability. In the 
present work, we focus on Monte Carlo simulations of 7, 
the second hyperpolarizability, to test the assumptions 
used in calculating the fundamental limits, to search for 
the broadest universal properties, and determine if this 
approach is capable of generating 7 values close to the 
limit. Finally, we will analyze the number of states that 
contribute to 7 when it approaches the limit - a first step 
in establishing the validity of the three-level ansatz for 
the second hyperpolarizability. 

II. APPROACH 

The second hyperpolarizability can be expressed using 
the dipole-free (DF) sum-over-states equation, which is 
equivalent to the traditional sum-over-states (SOS) ex- 
pression when enough states are included in the sum. 
The xxxx tensor component of the off resonant dipole- 
free SOS expression of the second hyperpolarizability is 
given by,[H 
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where a prime indicates that the ground state is excluded given by, 



in the summation. Eij = Ej 



Ej where E^ and Ej are 
Xij is the 



energies of i th and j th states, respectively, 
matrix element of position between states i and j. Note 
that Equation [2] holds for a three-dimensional molecule 
where we have chosen the coordinate system such that 
xxxx is the largest diagonal component of the second 
hyperpolarizability. 

The Thomas-Kuhn sum rules are the direct conse- 
quence of the commutator of the Hamiltonian with the 
position operator, [[H, x] ,x] where H is a Hamiltonian 
of the form 



V(t). 



(3) 



Here, p is the momentum at position r and V^(r) is the 
potential energy of the particle. The Thomas-Kuhn sum 
rules for the non-relativistic Hamiltonian Equation [3] arc 



E 

n=0 
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(4) 



where #,j is the Kronecker delta. Equation @] represents 
an infinite number of equations indexed by the integers 
{m,p). In the text that follows, we refer to a particular 
sum rule by stating the values (m,p). Note that Equation 
2] holds for a broader set of multi-electron Hamiltonians 
including spin, externally applied electromagnetic fields, 
and electron correlations. [23 Most of the cases can be 
described by the non-relativistic Hamiltonian. For the 
relativistic case, some modifications of Equation @] are 
required. [Hill 

The normalized form of the sum rules take the form 



E 



i 



(e, 



(5) 
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where ej = / E\q can take on values from one to infin- 
ity for i > 1 (the ground state is labeled by 0). Finally 
the normalized transition moments arc defined as 



where 
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Equations [7] and 2] imply that the largest possible tran- 
sition moment from the ground state is xqi because 
Eio < E i0 for i > 1. 

The three level ansatz is central to calculating the fun- 
damental limit of the second hyperpolarizability, 7-@ 
It states that when 7 is at the limit, only three states 
contribute. Using this ansatz and the sum rules given 
by Equation 31 the second hyperpolarizability given by 
Equation [2] is found to be constrained to the range, Q 
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Hence, defining the fundamental limit of the second hy- 
perpolarizability as 



= 4- 



~ 4 7i 4 
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then 



< lint < 1, 



(9) 



(10) 



where j int = y/~f max - Instead of studying 7 directly, 
we use the dimensionlcss intrinsic second hyperpolariz- 
ability because it takes into account simple scalings.[23j 
Hence, when 7 approaches e 4 N 2 h 4 /m 2 Ei , 7, n t ap- 
proaches unity. 

The numerical procedure is the same as in previous 
studies of the hyperpolarizability, 3:[Pj\ For a specific 
number of states s, the energy values are picked arbi- 
trarily and sorted in ascending order so that Eq < E\ < 
. . . < E s -i where Eq is the energy of the ground state 
and Ei is the energy of the i th excited state. For simplic- 
ity and with no loss in generality, we shift the energies so 
that the ground state energy is zero. Since = Eio/Eio, 
then eo = and e\ = 1. 

The next step is assigning transition moments that to- 
gether with the energies are forced to be consistent with 
sum rules. For a time-invariant Hamiltonian, it can be 
shown that the wavefunctions are real provided that there 
are no degeneracies. [27J Our method is not well suited to 
situations where degenerate states are present. As such, 
we treat the case where transition moments are real and 
intentionally avoid degeneracies. We are thus ignoring 
the class of problems in which (1) the Hamiltonian is not 
time invariant; and, (2) there are degeneracies. These 
cases will be considered in future studies. 



Since there are no dipole moments (the diagonal com- 
ponents of the transition moments) in jdf given by 
Equation^ we need only use the sum rules with m = p in 
Equation[5]to get all of the required transition moments. 
Starting with (m,p) = (0,0), we get 

eio|6o| 2 + e 2 o|6o| 2 + e 30 |6o| 2 + -.. + e„o|^o| 2 = 1. (11) 

Since e\o = ei — eo = 1, £01 can be randomly picked from 
the interval —1 < £01 < 1- Then £02 is obtained based 
on the fact that 



^ e no£ 2 o — 1 - e io£ 
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which implies 



£ 2 2 o < (1 - eiCi 2 ) /eao. 



(12) 



(13) 



A random number — 1 < r < 1 is used to get £20 from 
Equation [T3J 



60 = r^(l-ei^ )/e 20 . 



(14) 



For an s-state model the same method is used for all 
other transition moments, ^q, except for £, s -i,o, which is 
directly determined from the sum rules. For example, for 
a 10 state model, 



;io,o 



(1 



e 2 £ 2 



' e 9 £!o) /eio.o- (15) 



The remaining transition moments are calculated using 
the higher-order sum rules in sequence with p = 1,2,3,... 
in Equation[5j Since ^ is assumed to be real for all states 
i and j, 



(16) 



The energy and transition moment values are inserted in 
normalized form into Equation [2] to calculate 7^ . The 
normalized form is obtained by dividing Equation [5] by 
Equation [9] such that in Equation [21 7 — > 7j„t, Eij — > 
and Xij — > £ij . Note that one set of parameters obtained 
with the procedure given by Equations [T2] to Equation 
[TBI is called a single trial or single iteration. To iden- 
tify common properties of the second hyperpolarizability, 
hundreds of thousands of iterations are performed, and 
the distributions of the results are analyzed. 

Before proceeding, we need to address the issue of the 
signs of the transition moments. Recall that the diag- 
onal sum rules suffice to fix all the transition moments. 
However, all the transition moments in the diagonal sum 
rules appear as a squared modulus, . Thus, the di- 
agonal sum rules give no information about the sign of 
Xij. So one can choose the sign of each transition mo- 
ment independently from the values of all others. The 
question arises whether or not the signs are constrained 
by the sum rules that were neglected. 

To answer this question, we first refer to the non- 
diagonal sum rules for which m^pin Equation [SJ which 
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can be used to determine a condition on the signs by 
solving the equations for all different pairs of m and p. 
Considering the number of states, s, this approach de- 
mands solving s coupled equations. Since this is time 
consuming, we confronted the problem using the simpler 
approach of randomly picking all transition moments to 
be positive or negative. Neither /3 nor 7, calculated from 
the resulting matrix elements and energies, exceeded the 
fundamental limit. Thus we concluded that the sign of 
the transition moments can be picked randomly. We have 
studied the issue of truncated sum rules, [28[ and in fu- 
ture, we will address the self-consistency between the sum 
rules when truncated. 

The energies and transition moments are related via 
sum rules that impose a limit on all the transition mo- 
ments. To express the limit mathematically we start with 
Equation [5] with m = p, 

00 

$>„ p |£ np | 2 = 1. (17) 

For p = Equation [17] yields 
eio|6o| 2 + e 2 o|6o| 2 + e 3 o|C3o| 2 + - ■ - + e„o|e„ | 2 = 1. (18) 

Equation IT51 implies that the largest value of e^o l^io 1 2 
is attained when all other terms of the type e.,o|£jo| 2 (for 
j =/= i) are zero, i.e. 

e j0 |6o| 2 <l 0<i<n (19) 

or 

/1\V2 /1\ 1/2 

- — <6o< — • (20) 
\eioy \eioJ 

£io is largest when e^o is minimized and vice versa. Since 
based on our definition for energy spacing, when i > 1, 
we have e^o > 1 then £jo is limited to 

- 1 < 60 < 1 for < i < n. (21) 

Following the same procedure for p = 1 in Equation [3] 
results in 

eoi|Coi| 2 + e 2 i|6i| 2 + e 3 i|6i| 2 + -- - = 1, (22) 
whence 

e 2 i|6i| 2 + e 3 i|6i| 2 + --. + e„i|6 l i| 2 = 1 + eioieoil 2 - (23) 

Here we have used the fact that eoi = — eio. Substituting 
Equation 1191 into Equation 1231 gives an upper limit 

e 2 i|6i| 2 + e 3 i|6i| 2 + --. + e„i|£ni| 2 < 2. (24) 

This implies that eji|£ii| 2 < 2 (when 1 < i < n) and is 
obtained when for all j ^ i, eji\£ji\ 2 = 0. Then 

e*i|&i| 2 <2 Ki<n. (25) 



or 

/ 2 y/2 / 2 \ 1/2 

- — < 61 < — for 1< i < n. (26) 

Using the same method for p = 3,4, 5, . . . , n — 1, the 
constraint on the transition moments is found to be, 

2 p/2 2 P/ 2 

~^7J< k P for P < * < «• (27) 

p.- p.- 

Recall that we have assumed that the transition moments 
are real so £y = £ji . The sum rules specify the range of 
validity of all non-diagonal transition moments accord- 
ing to Equation [57] The Monte Carlo method typically 
yields energy spacings that lead to transition moments 
that lie below unity. Equation [57] is a general result that 
applied to quantum systems that are described by the 
Hamiltonian given by Equation [3] 



III. RESULTS AND DISCUSSION 

In this section, we use Monte Carlo simulations to 
study the statistics of the distribution of the intrinsic 
second hypcrpolarizability with the aim of understanding 
the physical properties of a quantum system that leads 
to a nonlinear response that approaches the fundamental 
limit. Correlations of the values of 7^4 with the average 
energy spacing and transition moments when 7i„ t is at 
the limit will be used to investigate whether such sys- 
tems behave in a way that suggests universal behavior. 
Finally, we discuss the effect that pairs of states have 
on the second hypcrpolarizability using a missing state 
analysis. 

When energies are chosen randomly, on average they 
are equally spaced. Thus the energy spacing is similar 
to the eigenenergies of a harmonic oscillator. To increase 
the domain of energy spacing that is probed by the Monte 
Carlo approach, we define an energy weighing factor /, 
such that 

Bk -> E{. (28) 

For / < 1 the energies are denser at higher energies and 
for / > 1 the states are denser at lower energies. 

Figure QJi shows the distribution of the second hyper- 
polarizability for a 15-state model for weight factors of 
/ = 0.75, 1, and 1.25. Frequency refers to the number of 
times 7i„t appears in the histogram. All simulated 7j„t 
values lie in the range predicted by Equation 1101 The 
maximum values of ji n t approaches unity, which suggests 
that the set of energies and transition moments that are 
consistent with sum rules can create an intrinsic second 
hyperpolarizabilities that can be close to fundamental 
limit (-fint — 1). On the other hand, the smallest j int 
peaks at zero and the largest negative value approaches 
—0.25, as predicted. 
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FIG. 1. (a)Distribution of simulated intrinsic second hyper- 
polarizabilities ji„t for a 15-state model and weight factors 
of 0.75, 1, 1.25. The inset shows histograms of 7;„ t near the 
limits. The bin size for the main figure is 0.003; and, for In- 
set b and Inset c are 0.05. The error bars are defined as the 
square root of the frequency. 



In contrast to studies that determine the shape of the 
class of potential energy functions that optimize the in- 
trinsic hypcrpolarizabilities, which leads to Pint < 0.708, 
the Monte Carlo method yields values of /3j n t that ap- 
proach the limit. We find that the same is true for ji n t- 
The distribution of ji n t for all three weighting factors 
peaks near zero. However, the distributions that are gen- 
erated with larger weighting factors have larger tails near 
the limits, as shown in the insets of Figure [TJ 

To better illustrate the distribution of 7 int near the 
limits, we have plotted the histogram of in the range 
[0.75, 1] (FigurcQJ)) and in the range [—0.1, —0.25] (Figure 
[Tb) where the bin size is 0.05. Both of these diagrams sug- 
gest that systems whose energy difference between adja- 
cent states gets larger for higher energy (black bars with 
/ = 1.25) appear more frequently near the limit. 

Recall that the three-level ansatz states that when the 
hyperpolarizability is near the fundamental limit, the 
ground and two excited states dominate. When the sum 
rules are applied to the three-level model of the hyper- 
polarizability, and the model is parameterized in terms 
of, 



_ £1 _ E 10 
E20 



< 2 



and 



where, 



X = & 
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(29) 
(30) 

(31) 



FIG. 2. a) Histogram of 7,:,^ for 5, 15 and 40 state models, 
(b) Enlarged histogram for ji n t £ [0.75, 1]. The bin size is 
0.05 



level ansatz and the fact that X = 1/ \/3 when the hyper- 
polarizability is large is found to be a universal property 
of all systems studied to date. [2^] While the fundamen- 
tal limit of the off-resonant second hyperpolarizability is 
calculated using the three-level ansatz, there is no ex- 
tensive body of literature that supports its validity. As 
we show below, the Monte Carlo simulations show that 
when "fi n t — > 1, the system can be modeled with three 
dominant states. 

Figure [2] shows the resulting distributions for Monte 
Carlo simulations of the intrinsic second hyperpolariz- 
ability for 5, 10, and 15 states. The inset (Figure [5b) 
focuses on the intrinsic hyperpolarizabilities in the range 
[0.75,1]. Clearly, when there are a fewer number of 
states, the frequency of ji nt near the limit increases sig- 
nificantly: the 5-state model yields an order of magni- 
tude higher frequency than the 15 and 40 state models. 
Also, the 15 state systems appear with higher frequency 
than the 40 state systems. This is a general trend that 
suggests the three level ansatz holds for the second hy- 
perpolarizability. 

Figure Q] shows that the largest second hyperpolariz- 
abilities are associated with greater energy spacing, sug- 
gesting that 



lim j int 



1. 



(32) 



the hyperpolarizability is then found to be at the funda- 
mental limit when X = l/</3and_E = 0.[l| The three- 



Figure [3] plots 7i„ t as a function of E for a 5 state model 
and one million iterations. When the intrinsic second 
hyperpolarizability approaches unity, the energy spacing 
most often falls in the range E < 0.1, suggesting that 
lint is the largest when E20 Eiq. the same result that 
is found for the first hyperpolarizability. 

To provide a more quantitative measure of the distri- 
bution of energies for various ranges of 7^4, we generate 
histograms from the data shown in Figure [3] The distri- 
bution of E for a 10 state model with 100,000 iterations is 
shown in Figured] Included are 12 equally spaced inter- 
vals of 7j„( (except for the range of (0.85, 1)). The points 
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FIG. 3. The distribution of 7^4 vs. E. For large values of 
the second hyperpolarizability, the majority of E values lie in 
the range [0, 0.1]. 
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FIG. 5. fi„t as a function of X for second hyperpolarizabilities 
within 5% of the limit. 
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FIG. 4. Distribution of E for different ranges of 7;„t for a ten- 
state model with 100,000 iterations. Each symbol corresponds 
to the specified range of jint- 



represent the Monte Carlo data and the curves are fits to 
a stretched exponential model of the form, 



In (F/F ) = - (E/Eo) 7 



(33) 



where Fq, Eq and n are fit parameters and F is the fre- 
quency. The stretched exponential was chosen as a model 
because it best fits the data with the fewest number of 
parameters. 

Note that all of the curves appear to be approxi- 
mately parallel to each other with the exception of the 
curves for j int in the range of (0.85, 1), (0.75,0.85), and 
(—0.25, —0.15); which fall off more steeply as a function 
of E. Thus, larger energy spacing is correlated with 7; nt 
near unity or near —0, 25, the positive and negative lim- 
its. 

Figure O shows 7, n t as a function of X for a 5-state 



FIG. 6. Distribution of X for different ranges of 7;, 



model and 100,000 runs when the second hyperpolariz- 
abilities are larger than 0.95. According to Figure [SJ as 
"tint approaches the limit, the range of X becomes nar- 
rower so that for the ji n t > 0.98, X lies in the range 
-0.2 < X < 0.2. 

To better quantify the results, Figure [6] shows the dis- 
tribution of X for various ranges of 7, n t using a 10-level 
model. The bin size is 0.05 and "/i n t has been divided into 
twelve equally-spaced ranges except for j int € [0.85,1]. 
lint € [—0.25,-0.15] is dominated by systems with 
X — ±1. For slightly smaller negative values of Jint, 
the peak in the distribution is less than unity and the 
rest of the curves peak at X = 0. 

While any value of 7^ is attainable for any arbi- 
trary choice of X, it is clear that the most likely value 
for the dominant transition moment is X = ±1 when 
lint — —0.25 and X — when ji nt = 1. Recall that 
Figure H] implies that the largest values of ji nt are for 
E = 0. Interestingly, the three-level ansatz shows that 
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when E = 0, X = yields ji n t = 1 and \X\ = 1 yields 
lint — — 0.25.0 Thus, while the the Monte Carlo simu- 
lations include systems where many states may be con- 
tributing to the second hyperpolarizability, on average, 
this statistical result gives the same result as the three- 
level ansatz. What remains to be investigated is the va- 
lidity of the the three- level ansatz. 

The three level ansatz is fundamental to the calcula- 
tions of the limit of the first- and second-order nonlinear 
optical response. The intrinsic hyperpolarizability /3 int , 
can expressed as fl7j . 

fknt = Yif%?, (34) 

n,m 

where the prime indicates that the ground state is ex- 
cluded in the summation. ^ s the fractional contri- 
bution of pairs of individual states n and m.[l7| 

f&? = (|) - ^F^) • (35) 

The total hyperpolarizability j3 can be calculated as 
the sum over all the fractional contributions. When the 
hyperpolarizability is near the fundamental limit, we find 
that /3™^™ is large for only one pair of states. This is in 
agreement with three level ansatz. 

Expressing "fi n t in terms of fractional contributions of 
pairs of states is complicated by the fact that three ex- 
cited states contribute to each term in Equation [2] We 
define l^ nt as the fractional contribution of the three 
states i, j and k 

lint = J]] J]] lint = J]] lint + 

i j k i^j^k i=j^k 

Y 7int+ Y lint+ Y ( 36 ) 
i^j=k i=k^j i=j=k 

In contrast to the hyperpolarizability, where /3™ l t = 0, 
the second hyperpolarizability does not vanish when i = 
j = k. 

The simplest approach for determining the contribu- 
tion of pairs of states is the missing state analysis, [29j 
where one calculates 7i„t in the absence of a particular set 
of states, {i, j), denoted by 7 l ™ ssm£ ' (i, j)- A comparison 
of lint with 7j™ ssms (i, j) describes the joint contribution 
of states i and j to the intrinsic second hyperpolarizabil- 
ity. The smaller the value of l™ t ssm9 (i, j), the larger the 
contribution of states i and j to ji n t ■ 

Figure [7] shows a logarithmic plot of the absolute value 
01 linV Sm9 {^ j) f° r a 10-state model with ji n t = 0.9886, 
the largest valued obtained in a Monte Carlo simula- 
tion. The most important pair of states are 1 and 2 with 
7 ™ ssm9 (2,l) = 2.76444 x 10~ 14 . The next most im- 
portant pair of states are 5 and 1 with 7j™ ssms (5, 1) = 

10- 11 . lTnt Sm9 (8, 1) and iZt Sm9 (9, 1) arc also on the or- 
der of 10 -11 . Thus, the dominant two states contribute 
about 10 3 times the next most important state. The 




FIG. 7. A log plot of |7int(i,i)| using the missing state anal- 
ysis for a 10 state model and yi„t = 0.9886. We remark that 
\lint(hj)\ = \lint{j,i)\ so that the two smallest values in the 
diagram are identical and represent only one pair of states. 

same procedure can be followed for other optimized val- 
ues, and similar results arc found. It is worth mentioning 

,1 , missinq / ■ -\ / missinq / ■ -\ 1 ,1 j- • j 

that 7 mt _ y {i,j) £ 7 int because the individ- 

ual terms in Equation [2] are not symmetric when two 
incidences are interchanged. Since each pair of states in 
Equation [5] contribute twice, we can define 

symmetrized/- -\ missinq / ■ -\ . missinq i ■ -\ fnn\ 

lint {h3)=lint (h])+lint {3,V, (37) 

for i ^= j and restrict the sum to i < j. We have applied 
such symmctrization to all values obtained with the miss- 
ing state analysis. 

IV. CONCLUSION 

Optimizing the intrinsic hyperpolarizabilities by vary- 
ing the potential energy function results in a 30% gap 
between the fundamental limit and optimized Pint- Ex- 
perimentally measured second hyperpolarizabilities also 
fall far-short of the limit for a wide range of molecules, 
but Monte Carlo simulations show that Pint can approach 
unity. In the present work, we find that Monte Carlo sim- 
ulations lead to values of 7 arbitrarily close to the funda- 
mental limit, with —0.25 < ji n t < 1, in agreement with 
analytical calculations. Q Also, the three level ansatz ap- 
pears to be obeyed as shown using the missing state anal- 
ysis, where the two dominant states account for over 99% 

Of lint- 

It is important to stress that the Monte Carlo approach 
may lead to a broader set of transition moments and en- 
ergies than are attainable with standard Hamiltonians; 
that is, Hamiltonians that include kinetic and potential 
energies of the many electrons in an atom as well as inter- 
actions with an external electromagnetic field. As such, 
universal behavior that may be typical of systems de- 
scribed by standard Hamiltonians may not be observed in 
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our calculations. Never-the-less, we find that the Monte 
Carlo simulations, on average, are consistent with ana- 
lytical results, with optimized -fi nt resulting when E « 
and \X\ = 1 for negative 7 and X = for positive 7. 

This is the first study to confirm that there may be 
universal properties associated with the second hyperpo- 
larizability when it is near the fundamental limit. As 



such, our approach may lead to the design of better ma- 
terials for third-order nonlinear-optical applications if the 
universal properties can be re-expressed in terms of pa- 
rameters that can be varied by a synthetic chemist. 
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